Introduction

This paper was developed out of the NSF-funded project Faunal Resource Depression and Intensification in the North American Southwest: Digital Data and Regional Synthesis (BCS-1153115). The project’s primary objective was to examine whether there was a relationship between the reduction in availability of large game and the intensification of turkey production in the late prehispanic period (1200–1500 CE) in the American Southwest. Addressing this topic required access to and the integration of multiple faunal data sets from Ancestral Pueblo sites dating to that time period. During the project these data sets were uploaded into tDAR (the Digital Archaeological Record, tdar.org), providing both the project personnel and the public access to the data. The faunal data sets curated during this project are available in the Archaeological Fauna: US Southwest collection on tDAR. We selected a subset of 32 of these data sets—those which had sufficient sample sizes for the taxa of interest—for this taphonomic analysis. The 32 datasets reported here are available in the Southwestern Taphonomic Protocol collection on tDAR.

Given that these assemblages had been analyzed and recorded by multiple faunal analysts with different coding schemes, an integrated analysis required that their variables be mapped to a set of general ontologies. This mapping was made possible by the development of such ontologies within tDAR. New software within tDAR then allowed the integrated analysis of these datasets (see Kintigh et al. (2018) for details on this component of the project).

Inter-site comparative analyses of the fauna data also necessitated an examination of the degree to which the zooarchaeological remains from these different sites had been affected by taphonomic processes (Bar-Oz and Dayan 2003; Gifford 1981; Lyman 1994; Lyman and Dean 2010). To assess variability in taphonomic history, Tiffany Clark (2014) developed a protocol that explores the degree to which different taphonomic factors played a role in assemblage formation. At the end of her portion of the project, we were concerned by the challenge of undertaking analyses of the diversity of taphonomic variables that were included in the protocol. Kyle Bocinsky thus joined the project to help develop this executable paper for the protocol so that it can be routinely used, or modified as necessary, by other researchers interested in the integrated analysis of faunal data. Clark (2014) contains details regarding the development of the protocol.

Overview

The remainder of this paper briefly discusses the development of the taphonomic protocol and executes it on the 32 faunal data sets included as part of our project. This analysis was developed as an R Markdown document — essentially, a text file that includes computer code in the R statistical language that runs all of the analyses discussed here and produces the figures and tables in this document. The HTML or PDF version of this paper (which you are probably reading) is “compiled” from that original text document and related data files, which were developed at https://github.com/bocinsky/swtp and are archived at https://doi.org/10.5281/zenodo.3256448. R Markdown allows for dynamic, data-driven analysis and report writing, and fosters reproducible research. As such, we hope that this document will not only serve as an introduction to the Southwestern Taphonomic Protocol, but will also be a template for others’ implementations of (or alterations to) it; we are currently developing an interactive tool to allow researchers to analyze their own data without having to change the R code. Details on using R Markdown can be found at http://rmarkdown.rstudio.com. The executable paper attempts to follow guidelines developed by Ben Marwick (2017) for practicing reproducible computing in archaeological research, and was developed using the rrtools package for R (Marwick 2019).

Overview of the SWTP

Selection of Variables and Quantification Method

The first step in constructing the taphonomic protocol was to identify the variables that could be used to assess the impact of different natural and cultural processes on assemblage formation. Data from the analysis of these variables could then be used to evaluate the relative degree of taphonomic comparability among assemblages and to identify individual data sets (or components therein) that display substantial taphonomic bias. A review of the literature found that a diversity of attributes has been employed to evaluate taphonomic biases in faunal assemblages (Bar-Oz and Dayan 2003; Bar-Oz and Munro 2004; Behrensmeyer, Allison and Briggs 1991; Lyman 1984, 1985, 1994; Marean 1991; Morlan 1994; Orton 2012; Pickering, Marean and Dominguez-Rodrigo 2003; Stiner 1992, 1994). Although much of this research has focused on examining a single or small set of variables (e.g., density-mediated attrition or fragmentation), several recent studies have sought to provide a more integrated approach to the study of faunal taphonomy (see Bar-Oz and Munro 2004; Marciniak 2001, 2005; and Orton 2012). These analytical schemes incorporate numerous variables in order to evaluate the effects that different taphonomic agents may have had on faunal assemblages.

The relatively comprehensive sets of variables that were employed by Bar-Oz and Munro (2004) provided a basis with which to begin to define the attributes that are most appropriate for the SWTP. However, those authors caution that the influence and combination of taphonomic factors that affect assemblages may differ significantly among time periods and across geographic regions, and as such it is important that researchers evaluate the suitability of individual taphonomic variables for their specific data sets and research questions. Towards this end, Clark (2014) identified a subset of the variables evaluated by Bar-Oz and Munro (2004) that appeared to be most applicable to the examination of taphonomic processes in the prehispanic Southwest. The variables that were selected could be evaluated using data that are commonly recorded on Southwestern sites by zooarchaeologists.

The variables fell into three broad analytical categories: bone surface modification, fragmentation intensity, and assemblage completeness. In this paper, these variables are compared across the three taxonomic categories that are relatively abundant in Southwestern assemblages and are of particular interest to the large game depletion/turkey intensification portion of the project: artiodactyls, lagomorphs, and turkey. Table 1 provides a list of the variables used in the protocol. Analyses assessing these variables have the potential not only to inform on the natural agents (weathering, gnawing, density-mediated attrition, and in situ attrition) that were involved in assemblage formation, but also on the degree of influence of anthropogenic factors involved in processing a carcass. A short description of each variable, including its interpretative potential and means of assessment, is provided in Table 1.

Table 1: List of possible taphonomic variables.
Variable Potential Taphonomic Information
Proportion of weathering of artiodactyl bone (NISP) Damage from natural peri-depositional formation processes
Proportion of gnawing on artiodactyl bone (NISP) Damage from natural peri-depositional formation processes
Proportion of burning of artiodactyl bone (NISP) Damage by human subsistence behaviors
Proportion of cut marks on artiodactyl bone (NISP) Damage by human subsistence behaviors
Degree of completeness of artiodactyl bone (NISP) Fragmentation from both natural and cultural agents
Artiodactyl average bone weight (grams) Fragmentation from both natural and cultural agents
Proportion of complete astragals (NISP) In situ attrition
Relationship between bone survivorship (NISP) and bone mineral density Density-mediated attrition
Frequency of burning of artiodactyl, lagomorph, and turkey bone (NISP) Intertaxonomic differences in damage by human subsistence behaviors
Frequency of cut marks on artiodactyl, lagomorph, and turkey bone (NISP) Intertaxonomic differences in damage by human subsistence behaviors
Degree of bone completeness for artiodactyls, lagomorphs, and turkeys (NISP) Intertaxonomic differences in natural and cultural agents
Bone survivorship (NISP) and bone mineral density of artiodactyls, lagomorphs, and turkeys Intertaxonomic differences in natural and cultural agents

During Clark’s development of the protocol (Clark 2014), although the Number of Identified Specimens (NISP) was used as the primary method of quantification, the Minimum Number of Individuals (MNI) was also calculated for select variables. Analyses of bone survivorship were undertaken using both NISP and MNI counts in order to evaluate the comparability of the results obtained with different quantification methods; both methods produced similar Spearman’s rank correlation coefficient values. These results suggest that in most cases, NISP counts can effectively be employed to assess the relationship between bone survivorship and bone mineral density. This conclusion supports the earlier findings by Grayson and Frey (2004), who determined that NISP-based body-part analyses could replicate the results of those obtained using derivative counting methods including MNI, MNE (Minimum Number of Elements), and MAU (Minimal Animal Unit). The NISP requires minimal time to calculate, can be easily replicated, and avoids the problems associated with the use of MNI counts on small-sized assemblages (see discussion in Grayson 2014). NISP is thus the quantification method used in the execution of the taphonomic protocol in this paper.

The final step in the development of the protocol was the delineation of the quantitative methods that may be used to gain a comparative understanding of how the taphonomic variables patterned across different sites. To this end, it was decided that in most cases, a visual assessment of the quantitative data was the most effective means for assessing taphonomic data across sites and for identifying those assemblages that may have been subject to taphonomic processes significantly different from the others. One of the strengths of this approach is that it allows the researcher to easily identify how, and to what extent, an assemblage may have been taphonomically altered. Based on this information, a researcher can determine whether or not an assemblage displays enough of a taphonomic bias to warrant its exclusion from a particular study. For the analysis of the degree of correlation between element representation and bone mineral density we used Spearman’s rank correlation coefficient, rho. For the rho data, we assessed both the numeric output and the graph. We bootstrap-estimated 95% confidence intervals around Spearman’s rho. In the graphs below, the dot denotes Spearman’s rho, while the line denotes the bootstrapped 95% confidence interval.

To visually compare the distribution of taphonomic data by variable, the proportion of each assemblage affected and the 95 percent confidence interval around the proportion are calculated in this paper. The dot denotes the proportion of each assemblage affected, with the black line indicating the error range at a 95 percent confidence level estimated using an exact binomial test (Clopper and Pearson 1934). The width of the confidence intervals can be used to assess the accuracy of the estimated proportion and the degree of confidence associated with the estimate. Narrow confidence intervals are often associated with large samples that display a low degree of variance; in contrast, wider confidence intervals tend to reflect more heterogeneous samples or samples that are relatively small in size. The use of the confidence intervals also allows the analyst to determine whether differences among assemblages are the result of taphonomic bias or are more likely attributable to sampling error.

Analysis and Discussion

Data overview

Data reported here are derived from publicly available faunal data sets from single-site contexts or from several sites within a local area. For the purpose of this paper, we chose not to assess taphonomic processes through time or to compare sites within data sets or proveniences within sites (though the SWTP can be used for both). Data sets were uploaded to tDAR either by the original faunal analysts or by project team members.

In several cases, including the Mesa Verde (Crow Canyon) and Salinas data sets, a single variable included multiple kinds of taphonomic data (e.g., burning, weathering, and gnawing). In order to include those sites in the taphonomic assessment we decided to disaggregate these data by creating additional variables in a copy of each data set. These copies are also available on tDAR but do not replace the original data sets there. The general caveat to this approach is that since only one of the several taphonomic conditions coded in a single variable could have been selected in the original analysis, elements that exhibited more than one condition would have had only one coded. In the case of Crow Canyon’s faunal database three modification columns were provided to the analyst although in most cases only one was used. This means that the results of the calculations undertaken in the executable paper are the minimum proportion of bones that were affected by a taphonomic condition for that site. We recommend future analyses code each type of taphonomic condition separately, preferably using a regionally-appropriate ontology such as the ontologies defined as part of the Archaeological Fauna: US Southwest collection on tDAR.

In a few other datasets a single taphonomic variable in our analysis was coded across multiple variables. For example, the Bryant Ranch faunal dataset includes separate columns for rodent gnawing and carnivore gnawing; we combined these into a new column mapped to the tDAR Faunal Gnawing Ontology. This was rare, but when it occurred we made a duplicate of the data set in tDAR and created a new synthetic variable that could be mapped to the ontology for that variable. Again, the original data set in tDAR remains intact.

The tDAR integration tool

The bases of tDAR’s integration tool are community-defined general ontologies to which variables from individual data sets are mapped (Kintigh et al. 2018). As part of the NSF project, members of the Southwestern Faunal Working Group defined 20 faunal ontologies including various taphonomic variables, which we used to integrate the 32 data sets in this analysis. We included ontologies for taxon, element, proximal/distal portion, burning intensity, butchering, completeness, gnawing, and weathering, plus a “count” column for recording NISP and a “display” column for recording weight.

Once the data sets were uploaded to tDAR and their relevant variables mapped (by the original analyst or project personnel) to the general ontologies in tDAR, we combined them using tDAR’s integration tool. The tDAR integration tool allows the export of the results of the integration for further analysis and to share with colleagues. Kintigh et al. (2018) provides a thorough overview of the tDAR integration tool; here, we detail how the integration for the SWTP was performed, including specific selections that must be made to reproduce our analysis.

In some data sets each line represented a single element; in others there was a “count” column that can be identified as such in tDAR. When a “count” column was present (Figure 1), we selected the variables that corresponded to NISP (as opposed to MNI). In most cases, this meant selecting “Actual #”, “N”, or something similar.

Screen capture of the tDAR integration tool count column selector. Users have to manually select the appropriate "count" variable, if one exists, to represent NISP.

Figure 1: Screen capture of the tDAR integration tool count column selector. Users have to manually select the appropriate “count” variable, if one exists, to represent NISP.

We then added a “display” column for weight (Figure 2), and selected the appropriate column for every data set in which specimen weight was recorded.

Screen capture of the tDAR integration tool display column selector, showing selected weight variables where available.

Figure 2: Screen capture of the tDAR integration tool display column selector, showing selected weight variables where available.

We then added integration columns (Figure 3). The tDAR integration tool allows the user to select categories within each ontology to be included in the integration, potentially collapsing lower-order categories into higher-order ones. For our analysis, however, we chose to handle the collapsing within R so as to have better control over collapsed categories. Therefore, when adding integration columns, users of the SWTP should click the “Select values that appear in any column” button in the left hand column. This will select the check boxes for categories represented in any of the integrated datasets.

Screen capture of the tDAR integration tool integration column selector. For integration columns in our analysis, users should click the "Select values that appear in **any column**" button in the left hand column.

Figure 3: Screen capture of the tDAR integration tool integration column selector. For integration columns in our analysis, users should click the “Select values that appear in any column” button in the left hand column.

Once all integration columns were selected (Figure 4), we requested the tDAR servers to process the integration. Due to the large amount of data in these datasets, the integration tool timed out and we had to request the integration output from the tDAR staff directly (they happily provided it many, many times). The raw output from our tDAR integration reported in this paper is available here. The processed data are reported in Supplementary Table 1. All tabular data presented here and in supplementary information can be downloaded here.

Screen capture of the tDAR integration tool showing the count, display, and integration columns.

Figure 4: Screen capture of the tDAR integration tool showing the count, display, and integration columns.

Collapsing taxon data into artiodactyls, lagomorphs, and turkeys

The SWTP focuses on three taxonomic categories: artiodactyls, lagomorphs, and turkeys. Given that many datasets had identified taxa to more specific levels than order, we developed two new packages for R. The first, called tdar, allows R users to access the tDAR Appplication Programming Interface (API), enabling credentialed users to download datasets on tDAR by their tDAR dataset identifier. Here, we use the tDAR package to automatically download the ontologies in the Archaeological Fauna: US Southwest collection on tDAR, though other users might find the package particularly useful for downloading archived archaeological datasets. Access to the tDAR API is currently by request only; please refer to the README document for the tdar package for information on gaining access.

The second R package, treecats, manipulates tree-like categorical data such as the tDAR ontologies and enables users to collapse nested categories in a comprehensive and reproducible way. While the tDAR integration tool allows for collapsing ontologies, it currently depends on the user selecting the appropriate higher-order check boxes, a process that in our experience was tedious and perhaps too error prone to be sufficiently reproducible. It is possible, however, to save a specific integration within tDAR to run again on the same or different data sets. Here the treecats package parses the Archaeological Fauna: US Southwest ontologies and allows for on-the-fly ontological filtering and collapsing.

Using treecats, we collapsed all species in the taxonomic order Artiodactyla, all the species in the taxonomic order Lagomorpha, and included any avian remains categorized as “large aves” in Meleagris gallopavo. Our goal was to conform to established precedent in southwestern archaeology, while using sufficiently broad taxonomic classes so as to account for analyst error. Although it is somewhat common in the US Southwest to include “large mammal” in analyses involving artiodactyl fauna, we did not include large mammal counts in this analysis because some datasets in our analysis only included identified fauna and not elements from this more general category. Thus the results of our calculations would have been incommensurate across sites if large mammal counts were included. For the bone survivorship analysis, we also collapsed the proximal/distal portion data into “Proximal”, “Distal”, or “Complete”.

Data summary

Table 2 summarizes the data sets, including the counts of samples with valid measurements across each variable reviewed in the Southwest Taphonomic Protocol. Our final dataset includes a total NISP of 329,379 across all taxa, including 41,412 artiodactyl specimens, 98,127 lagomorph specimens, and 29,336 turkey or large bird specimens.

Table 2: Summary of datasets and variables included in this analysis.

Thresholds for dataset inclusion

In the initial stages of developing this executable paper we found that there were assemblages in which relevant variables were present, but rarely used, or where sample sizes were far too small to result in meaningful information. Particularly in the case of the rare use of a variable, calculations of proportions and confidence intervals resulted in potentially spurious output. We thus established thresholds for inclusion of variables in the taphonomic protocol that would preclude their calculation where results would be misleading. The threshold for inclusion of artiodactyls, lagomorphs, or turkey elements in the analysis of a variable was set at an NISP of 25 for the taxon. Additionally, 40 percent of the taxon assemblage had to have been coded for that variable. Thus, the minimum NISP used in the calculation of the proportion of a taxon assemblage affected by a particular taphonomic variable is ten (40% of 25). When assessing assemblage completeness, we only included datasets with a minimum of 5 NISP for the element under consideration (astragali and 3rd phalanxes). Our bone surviorship analysis includes all datasets with a minimum of 25 total NISP among the elements under consideration.

Related to thresholds, a number of assemblages used the coding option of “indeterminate.” While this meant that the variable was indeed coded for all of the cases, the actual content of that coding did not contribute information to the interpretation of how prevalent the taphonomic condition was in that assemblage. We thus chose to consider those cases uncoded for analyzing taphonomic differences across sites. Therefore sites with more than 60 percent of the cases for each taxon coded as “indeterminate” fell below the minimum of threshold of 40 percent coded and were not included in the calculation of that variable.

Table 3 shows which datasets were included in each analysis for each taxon after thresholding. Table provides the number of datasets included in each analysis. Note that the assemblage completeness and bone survivorship analyses are not included in these tables.

Table 3: Datasets and variables included in this analysis after thresholding.

Bone surface modification

Burning

Proportion of assemblage showing evidence of burning
Burning is often a result of food preparation or refuse disposal practices and as such, can be an important cultural agent in assemblage formation. However, because burning renders bones more susceptible to fragmentation (Stiner et al. 1995), this variable may also provide insight into post-depositional processes. Burning is assessed by calculating the proportion of specimens (NISP) in the artiodactyl, lagomorph, and Meleagris gallopavo assemblages showing signs of exposure to heat (e.g., partial charring, charring, or calcined); these were mapped to “Burned” or “Probably Burned” in the SWUS Fauna Burning Intensity Ontology.

Select a taxon to view its burning data:

Artiodactyla

Figure 5: Artiodactyl burning. Data provided in Supplemental Table 3.

Lagomorpha

Figure 6: Lagomorph burning. Data provided in Supplemental Table 4.

Meleagris gallopavo

Figure 7: Meleagris gallopavo burning. Data provided in Supplemental Table 5.

Butchering

Proportion of assemblage showing butchering marks
The occurrence cut marks in a faunal assemblage may directly inform on butchering practices and processing techniques. The analysis of cut marks is used in the taphonomic study to evaluate the extent of damage resulting from human subsistence behavior. This variable is assessed by calculating the proportion of specimens (NISP) in the artiodactyl, lagomorph, and Meleagris gallopavo assemblages that show evidence of butchering or cut marks; these were mapped to “Butchered” or “Probably Butchered” in the SWUS Fauna Butchering Ontology.

Select a taxon to view its butchering data:

Artiodactyla

Figure 8: Artiodactyl butchering. Data provided in Supplemental Table 6.

Lagomorpha

Figure 9: Lagomorph butchering. Data provided in Supplemental Table 7.

Meleagris gallopavo

Figure 10: Meleagris gallopavo butchering. Data provided in Supplemental Table 8.

Gnawing

Proportion of assemblage showing evidence of gnawing
The extent of rodent and carnivore gnawing in an assemblage may provide direct information regarding peri-depositional formation processes. This variable is assessed by calculating the proportion of specimens (NISP) in the artiodactyl, lagomorph, and Meleagris gallopavo assemblages that were mapped to “Gnawed” in the SWUS Fauna Gnawing Ontology.

Select a taxon to view its gnawing data:

Artiodactyla

Figure 11: Artiodactyl gnawing. Data provided in Supplemental Table 9.

Lagomorpha

Figure 12: Lagomorph gnawing. Data provided in Supplemental Table 10.

Meleagris gallopavo

Figure 13: Meleagris gallopavo gnawing. Data provided in Supplemental Table 11.

Weathering

Proportion of assemblage showing evidence of weathering
The extent of weathering within an assemblage may provide direct information regarding peri-depositional formation processes. Weathering is assessed by calculating the proportion of specimens (NISP) in the artiodactyl, lagomorph, and Meleagris gallopavo assemblages that were mapped to “Weathered” in the SWUS Fauna Weathering Ontology.

Select a taxon to view its weathering data:

Artiodactyla

Figure 14: Artiodactyl weathering. Data provided in Supplemental Table 12.

Lagomorpha

Figure 15: Lagomorph weathering. Data provided in Supplemental Table 13.

Meleagris gallopavo

Figure 16: Meleagris gallopavo weathering. Data provided in Supplemental Table 14.

Fragmentation intensity

Fragmentation

Proportion of assemblage that is highly fragmented
The proportion of highly fragmented bone can be used to evaluate the overall degree of fragmentation among assemblages. As noted by Todd and Rapson (1988: 33–35), a variety of different processes may result in greater fragmentation of bone elements within an assemblage. These may include pre-depositional factors related to grease and marrow extraction or post-depositional damage resulting from overburden/sedimentary compaction or trampling. The degree of fragmentation of artiodactyl, lagomorph, and Meleagris gallopavo bone was assessed in this study by calculating the proportion of bone within each assemblage that was less than 25 percent complete; these were mapped to “<25%” in the SWUS Fauna Completeness Ontology.

Select a taxon to view its fragmentation data:

Artiodactyla

Figure 17: Artiodactyl fragmentation. Data provided in Supplemental Table 15.

Lagomorpha

Figure 18: Lagomorph fragmentation. Data provided in Supplemental Table 16.

Meleagris gallopavo

Figure 19: Meleagris gallopavo fragmentation. Data provided in Supplemental Table 17.

Weight

Bone weight
There is a broad correlation between fragment size, count, and weight (Lyman 2008: 102–103). Calculating the median bone weight for taxa can provide a means for evaluating the degree of fragmentation associated with a specific assemblage. Here, we show the mean, median, and interquartile range (IQR) of bone weights. In the graph below, the closed circles represent the medians, and the open circles represent the means. Unsurprisingly, the weight distributions are universally right skewed, with low numbers of heavier elements increasing the mean well above the median bone weight.

Artiodactyla

Figure 20: Artiodactyl bone weight. Data provided in Supplemental Table 18.

Lagomorpha

Figure 21: Lagomorph bone weight. Data provided in Supplemental Table 19.

Meleagris gallopavo

Figure 22: Meleagris gallopavo bone weight. Data provided in Supplemental Table 20.

Assemblage completeness

Astragal completeness

Proportion of complete artiodactyl astragals
Marean (1991) has argued that the extent of completeness of certain dense bone elements, such as carpals and tarsals, within an assemblage may be used to indirectly measure the degree of in situ attrition. Specifically, he posits that because these elements have little or no grease or marrow, they should not be subject to bone fragmenting behaviors by humans or carnivores. Rather, a high degree of fragmentation of large mammal carpals and tarsals more likely reflects the post-burial destruction of bones. Marean (1991: 687) further notes that if elements as dense as the astragalus or calcaneus are severely damaged within an assemblage, then it is likely that less dense bones have been completely destroyed. To evaluate in situ attritional processes in the current study, the proportion of complete artiodactyl astragals was calculated for each assemblage using the SWUS Fauna Completeness Ontology as above.

Figure 23: Artiodactyl astragals. Data provided in Supplemental Table 21.

3rd Phalanx completeness

Proportion of complete artiodactyl 3rd phalanxes
Following Marean (1991) logic, noted for the astragali, a second completeness variable was included in the protocol—the proportion of completeness for the slightly less dense third phalanx [0.25 bone density value for Odocoileus spp.]. This addition allows for more detailed examination of the taphonomic patterns associated with in situ bone attrition. The proportion of complete 3rd phalanxes was calculated for each assemblage using the SWUS Fauna Completeness Ontology as above.

Figure 24: Artiodactyl phalanxes. Data provided in Supplemental Table 22.

Bone survivorship and mineral density

Artiodactyl bone survivorship and bone mineral density
A number of taphonomic studies have shown that bone mineral density may be an important factor in skeletal element representation (e.g., Binford 1981; Kreutzer 1992; Lyman 1984, 1994; 1992). This research indicates that certain parts of a skeleton are consistently better preserved in archaeological assemblages due to their high structural density (Lyman 1994: 234–258). To determine if elemental representation within any of the faunal assemblages included in this study had been affected by differential density-mediated attrition, NISP values for 16 different artiodactyl element bone parts—femur (distal/proximal), humerus (distal/proximal), radius (distal/proximal), tibia (distal/proximal), ulna (distal/proximal), pelvis, astragalus, atlas, calcaneus, mandible, and scapula—were compared with bone structural density data for Odocoileus spp. (Lyman 1984, Table 6, value \(VD\)). The strength of the relationship between element abundance and bone density was assessed using Spearman’s rank correlation coefficient.

Figure 25: Artiodactyl bone survivorship. Data provided in Supplemental Table 23.

Intertaxonomic comparison

Burning

Proportion of burning of artiodactyl, lagomorph, and turkey bone
This variable compares the pattern of burning among subgroups in the assemblage in order to determine if artiodactyl, lagomorph, and turkey remains were differentially affected by food preparation or refuse disposal behaviors. Burning is assessed by calculating the proportion of artiodactyl, lagomorph, and turkey specimens (NISP) in the assemblages that show evidence of exposure to heat (e.g., partial charring, charring, or calcined).

Figure 26: Intertaxonomic burning. Data provided in Supplemental Table 24.

Butchering

Proportion of butchering marks on artiodactyl, lagomorph, and turkey bone
This variable compares the proportion of butchering marks among artiodactyl, lagomorph, and turkey remains to examine intertaxonomic variability in the extent of damage resulting from human subsistence behavior. This variable is evaluated by calculating the proportion of artiodactyl, lagomorph, and turkey specimens (NISP) in each artiodactyl that shows evidence of butchering or cut marks.

Figure 27: Intertaxonomic butchering. Data provided in Supplemental Table 25.

Fragmentation intensity

Degree of fragmentation of artiodactyls, lagomorphs, and turkeys
To assess intertaxonomic variability in fragmentation patterns, the proportion of complete bone within the assemblage was calculated for artiodactyls, lagomorphs, and turkeys. Significant variation in proportion values among these groups may indicate that taphonomic agents differently affected the bones of specific taxon. This variable is determined by calculating the proportion of bone that more than 50 percent complete (NISP) within the artiodactyl, lagomorph, and turkey assemblages.

Figure 28: Intertaxonomic fragmentation. Data provided in Supplemental Table 26.

Bone survivorship

Bone survivorship and bone mineral density of artiodactyls, lagomorphs, and turkeys
This variable compares the relationship between bone survivorship and bone mineral density among taxonomic subgroups in order to assess the differential effects of density-mediated attrition on artiodactyl, lagomorph, and turkey remains. To determine if certain subgroups were differentially affected by this bias, NISP values for a number of bone element parts were calculated for artiodactyl, lagomorph, and turkey. The abundance of these elements was then compared with bone structural density data from Odocoileus spp. (Lyman 1984, Table 6, value \(VD\)), leporid (Pavao and Stahl 1999, Table 1, value \(VD_{LD/BT}\)), and Meleagris gallopavo (Dirrigl 2001, Table 1, value \(BMDa\))), respectively. The strength of the relationship between element abundance and bone density for each subgroup was assessed using Spearman’s rank correlation coefficient.

Figure 29: Intertaxonomic bone survivorship. Data provided in Supplemental Table 27.

Discussion

In this paper we immersed ourselves in the application of the Southwestern Taphonomic Protocol (SWTP) to thirty-two late prehispanic faunal assemblages. Our goal has been to compare the taphonomic histories of these datasets to determine whether any are too dissimilar to include in subsequent substantive analyses of anthropological questions. Application of the SWTP has allowed us to identify several broad patterns that we consider here.

Variable Use and Utility

In developing the original Southwestern Taphonomic Protocol, Clark (2014) drew from global faunal literature and several pilot analyses to identify those variables most likely to be effective in distinguishing different aspects of taphonomic histories in southwestern faunal assemblages. In our implementation we found that some of these variables are more useful simply because they tend to be coded by more analysts. More engagement with taphonomic analyses might encourage zooarchaeologists to expand their coding to include at least some of the other variables used in the protocol. In addition, as we noted earlier in this paper, in some datasets relevant taphonomic data had been coded in a catchall variable that included multiple factors such as butchering, gnawing, and weathering. Coding these variables separately would make them more amenable to taphonomic analysis. The fact that they are already analyzed in these datasets should make it relatively straightforward to disaggregate the catchall variable since variable states have already been identified.

We have realized that, given the relative lack of attention to taphonomic variables, at least in the US Southwest, it is not surprising that bone mineral density (used in our bone survivorship analysis) has emerged over the decades as a dominant variable in taphonomic analyses. The basic data (element) that this variable requires are available in most datasets. Table 2 presents the number of datasets for which there were different types of taphonomic data.

The least available SWTP data involved specific elements, artiodactyl astragali and 3rd phalanx. Sample sizes were often quite small and thus the confidence intervals are generally too large for their fragmentation rates to be useful in distinguishing among sites.

Sample size is also an issue for the bone mineral density and bone survivorship computations. Confidence interval ranges are very broad for this comparative analysis because it requires counts of specific identifiable elements, which especially for artiodactyls are few in number. Most artiodactyl bone in southwestern archaeological sites is shaft fragments. Moreover, because the analysis requires proximal and distal portion information, it was only possible to undertake it for twenty of the thirty-two assemblages. The small sample sizes associated with specific artiodactyl elements likely influenced the p values generated in the bone mineral analysis. The p values are all non-significant except for Sand Canyon and Yellow Jacket, for which the relationship between bone mineral density and bone survivorship is highly significant.

Variable-wise Distinctions among Sites

Burning: For the most part, southwestern faunal assemblages exhibit relatively low frequencies of burning, generally below 20% for all taxonomic categories we examined. The one site in these assemblages that stands out for relatively high burning of artiodactyl and lagomorph bone is Frank’s Ruin, which itself was attacked and burned (Solometo, Rautman and Chamberlin 2017). Overall, the highest burning rates are for artiodactyl bone in general, with lagomorph and turkey less so. The distribution of artiodactyl burned frequencies makes the point that it is the analyst who will determine which sites are too taphonomically different to include in an analysis. Is a 15% burning threshold more appropriate (excluding sites from Albert Porter and above in terms of artiodactyl burning frequency), or should the threshold be above 30% burning (excluding Alameda School and above)?

Butchering: Southwestern fauna rarely show evidence of butchering. In all thirty-two datasets it appears in less than 10% of faunal elements and thus is not a distinguishing variable for the late prehispanic US Southwest.

Gnawing: Gnawing was also recorded on less than 10% of faunal elements, with the exception of the slightly higher frequencies on Sand Canyon and San Antonio de Padua. Thus, it too does not appear to be a distinguishing variable for taphonomic histories in the late prehispanic US Southwest.

Weathering: Weathering of faunal bone is also relatively uncommon, at less than 10% frequency for the three taxonomic categories in most of these assemblages. This is not surprising given the dry southwestern climate. Bryant Ranch and Pena Blanca, however, are clearly outliers with weathering being particularly high for lagomorphs and artiodactyls, and Pena Blanca for turkey as well. San Antonio exhibits a relatively high frequency of weathered artiodactyl bone. The high frequency of weathered turkey bone at Gran Quivira is at odds with the other data and may be an artifact of the context in which they were recovered—largely burials as opposed to midden. As we discuss below, high weathering frequencies for Pena Blanca and San Antonio de Padua are likely an artifact of the two-variable coding scheme used in recording the fauna from this site. Note that we have not calculated weathering for lagomorphs at Bryant Ranch as the sample size is below 25.

Fragmentation intensity: The distribution of fragmentation rates is relatively continuous from 30-50%, again making it up to the analyst to decide what the threshold should be for including or excluding faunal assemblages in an analysis based on their degree of fragmentation. A few sites, however, do stand out in terms of higher fragmentation rates Pena Blanca and San Antonio assemblages exhibit high frequencies of fragmented artiodactyl bone, and the rates for HARP and Kite somewhat lower. Pena Blanca lagomorph bone is also highly fragmented, as is San Antonio turkey bone. Artiodactyl bone is more fragmented than lagomorph bone in general, which is not surprising given the more intense butchering of larger game.

Astragal completeness: As noted above, the wide confidence intervals (due to low sample sizes) make this variable are not terribly helpful, but the data do reinforce the relatively high fragmentation of the Pena Blanca artiodactyl assemblage. No other fragmentation analyses, however, indicate that the Pueblo Colorado faunal assemblage is more fragmented than those from the other sites.

3rd Phalanx completeness: Again, the confidence intervals are wide, but San Antonio does appear to be more fragmented.

Artiodactyl bone survivorship: As noted above, only the results for Sand Canyon and Yellow Jacket yielded significant results when calculating Spearman’s rank correlation between bone survivorship and bone mineral density. It is unlikely that artiodactyl bone survivorship is useful without substantially larger sample sizes of artiodactyl elements.

Interpretation of variable output

The immense value of the diversity of variables in the SWTP is that site assemblages have been taphonomically altered by different variables. Thus focusing on only one or two would miss important taphonomic processes. The data also allow analysts to draw more fine-grained distinctions in taphonomic alteration among sites should they wish to do so by selecting specific thresholds in the frequency of different kinds of alteration beyond which they believe assemblages are too modified, or by investigating whether a particular variable has a significant impact on the interpretative value of an assemblage.

Table 4: Sites identified as outliers in the SWTP analyses. Bolded assemblages are the primary outliers.
Variable Sites
Burning Franks Ruin, RCAP, Homol’ovi IV, Pena Blanca, Alameda School, LA 12587 (Lagomorpha)
Butchering none
Gnawing Sand Canyon Pueblo, San Antonio de Padua (LA 24)
Weathering Bryant Ranch, Pena Blanca, San Antonio
Fragmentation Pena Blanca, San Antonio, LA 12587 (Lagomorpha)
Bone survivorship Yellow Jacket Pueblo, Sand Canyon Pueblo

Two of the strong outliers in the weathering portion of the analysis, however, alert us to the importance of understanding the nature of the coding scheme the analyst used before assuming that marked differences among sites are behaviorally meaningful. In the case of Pena Blanca and San Antonio de Padua, weathering was coded across two variables, ‘environmental alteration’ which identified the type of weathering on the bone (root etching, polished, corroded, etc.) and ‘environmental degree,’ which was coded from light to heavy. At both sites there were high frequencies of ‘root etching’ that was ‘light.’ In consulting with the original faunal analyst, Nancy Akins, she is of the opinion that this kind of weathering was likely not recorded as weathering by other faunal analysts. The apparent high frequency of weathering recorded at those two sites, then, is probably the result of the more fine-grained recording of weathering than actual differences in the proportion of weathered bone from those recorded at most other sites.

The analysis of these multiple variables also suggests, however, that some that were in the original taphonomic protocol can be removed for late period southwestern faunal assemblages. Astragal and 3rd phalanx completeness information is redundant with the results of the broader fragmentation analysis, and given the sample size issue with these two variables it does not appear to be necessary to include them. Similarly, the relative lack of evidence for butchering and for gnawing on bone, at least for these kinds of southwestern assemblages, indicates that these two variables would not need to be included either.

Intertaxonomic variation

The intertaxonomic graphs allow for a visual inspection of variation in taphonomic history among taxonomic categories both within and among sites. The burning, butchering, and bone survivorship charts provide a good summary that visually brings together the data considered separately in the previous section of the paper. The graphs allow for the more in-depth study of the differential impact that taphonomic factors may have on the taxonomic groups within an assemblage. For example, an examination of the Burning graph indicates that several assemblages, including the Alameda School Site, Franks Ruin, Homolovi IV, and Pena Blanca, all exhibit much higher frequencies of burning among artiodactyl bone compared to either lagomorph or turkey remains. Based on these findings, it appears that burned artiodactyl remains are largely responsible for the elevated level of burning observed in these assemblages.

The Fragmentation comparative chart is particularly important for thinking about substantive comparative analyses among these faunal assemblages. That graph makes clear that the three taxonomic categories vary in fragmentation rates within sites and that this variation is by no means constant across sites. It is thus possible that differential fragmentation rates within and across sites will affect the output from the calculation of ratio indices, such as the artiodactyl index.

Conclusions

In conclusion, the application of the Southwestern Taphonomic Protocol to over thirty late prehispanic southwestern faunal datasets has allowed us to identify those sites whose taphonomic histories appear considerably different from the majority of datasets in our study. We hope that these important results will spur fellow zooarchaeologists to include a diversity of taphonomic variables in their data collection both to inform their own analyses of faunal utilization in their areas, and to expand our current understandings of taphonomic processes in the US Southwest.

References

Bar-Oz, G and Dayan, T 2003 Testing the use of multivariate inter-site taphonomic comparisons: The faunal analysis of Hefzibah in its Epipalaeolithic cultural context, Journal of Archaeological Science, 30(7): 885–900.
Bar-Oz, G and Munro, N D 2004 Beyond cautionary tales: A multivariate taphonomic approach for resolving equifinality in zooarchaeological studies, Journal of Taphonomy, 2(1): 201–221.
Behrensmeyer, A K, Allison, P A and Briggs, D E 1991 Terrestrial vertebrate accumulations. In: Allison, P.A. and Briggs, D.E. (eds.) Taphonomy: Releasing the data locked in the fossil record. London: Springer. pp. 291–335.
Binford, L R 1981 Bones: Ancient men and modern myths. San Diego, California: Academic Press.
Clark, T C 2014 A Faunal Taphonomic Protocol for the Southwestern US.
Clopper, C J and Pearson, E S 1934 The use of confidence or fiducial limits illustrated in the case of the binomial, Biometrika, 26(4): 404–413.
Dirrigl, F J 2001 Bone mineral density of wild turkey (meleagris gallopavo) skeletal elements and its effect on differential survivorship, Journal of Archaeological Science, 28(8): 817–832.
Gifford, D P 1981 Taphonomy and paleoecology: A critical review of archaeology’s sister disciplines, Advances in Archaeological Method and Theory, 4: 365–438.
Grayson, D K 2014 Quantitative zooarchaeology: Topics in the analysis of archaelogical faunas. New York: Academic Press.
Grayson, D K and Frey, C J 2004 Measuring skeletal part representation in archaeological faunas, Journal of Taphonomy, 2(1): 27–42.
Kintigh, K W, Spielmann, K A, Brin, A, Candan, K S, Clark, T C and Peeples, M 2018 Data integration in the service of synthetic research, Advances in Archaeological Practice, 6(1): 30–41.
Kreutzer, L A 1992 Bison and deer bone mineral densities: Comparisons and implications for the interpretation of archaeological faunas, Journal of Archaeological Science, 19(3): 271–294.
Lyman, R L 1984 Bone density and differential survivorship of fossil classes, Journal of Anthropological Archaeology, 3(4): 259–299.
Lyman, R L 1985 Bone frequencies: Differential transport, in situ destruction, and the MGUI, Journal of Archaeological Science, 12(3): 221–236.
Lyman, R L 1994 Vertebrate taphonomy. Cambridge, England: Cambridge University Press.
Lyman, R L 2008 Quantitative paleozoology. Cambridge, England: Cambridge University Press.
Lyman, R L and Dean, R M 2010 Prehistoric anthropogenic impacts to local and regional faunas are not ubiquitous. In: Dean, R. (ed.) The archaeology of anthropogenic environments. Center for archaeological investigations occasional papers. Carbondale, Illinois: Southern Illinois University. pp. 204–224.
Lyman, R L, Houghton, L E and Chambers, A L 1992 The effect of structural density on marmot skeletal part representation in archaeological sites, Journal of Archaeological Science, 19(5): 557–573.
Marciniak, A 2001 Scientific and interpretive components in social zooarchaeology: The case of early farming communities in Kujavia, Archaeologia Polona, 39: 87–110.
Marciniak, A 2005 Social changes in the early European Neolithic: A taphonomy perspective. In: O’Connor, T.P. (ed.) Biosphere to lithosphere: New studies in vertebrate taphonomy. Oxford, England: Oxbow. pp. 148–154.
Marean, C W 1991 Measuring the post-depositional destruction of bone in archaeological assemblages, Journal of Archaeological Science, 18(6): 677–694.
Marwick, B 2017 Computational reproducibility in archaeological research: Basic principles and a case study of their implementation, Journal of Archaeological Method and Theory, 24(2): 424–450.
Marwick, B 2019 rrtools: Creates a Reproducible Research Compendium.
Morlan, R E 1994 Bison bone fragmentation and survivorship: A comparative method, Journal of Archaeological Science, 21(6): 797–807.
Orton, D C 2012 Taphonomy and interpretation: An analytical framework for social zooarchaeology, International Journal of Osteoarchaeology, 22(3): 320–337.
Pavao, B and Stahl, P W 1999 Structural density assays of leporid skeletal elements with implications for taphonomic, actualistic and archaeological research, Journal of Archaeological Science, 26(1): 53–66.
Pickering, T R, Marean, C W and Dominguez-Rodrigo, M 2003 Importance of limb bone shaft fragments in zooarchaeology, Journal of Archaeological Science, 30(11): 1469–1482.
Solometo, J, Rautman, A and Chamberlin, M 2017 Social risk and conflict among the plaza-pueblo villages of the Salinas province 1100–1400 CE. In: Spielmann, K.A. (ed.) Landscapes of social transformation in the salinas province and the eastern pueblo world. Tucson, Arizona: University of Arizona Press. pp. 69–102.
Stiner, M C 1992 Overlapping species ‘choice’ by Italian Upper Pleistocene predators, Current Anthropology, 33(4): 433–451.
Stiner, M C 1994 Honor among thieves: A zooarchaeological study of Neandertal ecology. Princeton, New Jersey: Princeton University Press.
Stiner, M C, Kuhn, S L, Weiner, S and Bar-Yosef, O 1995 Differential burning, recrystallization, and fragmentation of archaeological bone, Journal of Archaeological Science, 22(2): 223–237.
Todd, L C and Rapson, D J 1988 Long bone fragmentation and interpretation of faunal assemblages: Approaches to comparative analysis, Journal of Archaeological Science, 15(3): 307–325.